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MULTI-GRID  CALCULATION  OF  THREE-DIMENSIONAL 


TRANSONIC  POTENTIAL  FLOWS 


D.  A.  Caughey* 

Cornell  University 
Ithaca,  New  York  14853 


ABSTRACT 

A  multi-grid  algorithm  has  been  developed  to 
speed  the  iterative  convergence  of  calculations 
for  the  transonic  potential  flow  past  swept  wings 
and  wing-fuselage  combinations.  The  method  is 
based  upon  a  fully-conservative,  finite-volume 
approximation  to  the  steady  potential  equation 
which  is  second-order  accurate  everywhere  in  the 
flow  field  except  near  shock  waves.  The  multi¬ 
grid  scheme  is  incorporated  within  the  framework 
of  an  alternating  successive-line-overrelaxation 
(SLOR)  solver  of  the  difference  equations.  Com¬ 
puted  results  confirm  the  second-order  accuracy 
of  the  scheme,  and  demonstrate  the  effectiveness 
of  the  multi-grid  procedure. 


I.  INTRODUCTION 

In  the  past  several  years,  algorithms  have 
been  developed  for  predicting  the  transonic 
potential  flow  past  reasonably  complete  aircraft 
configurations.  In  particular,  the  finite-volume 
awthod  of  Jameson  and  Caughey 1-3  has  made  it  pos¬ 
sible  to  calculate  the  transonic  potential  flow 
past  any  configuration  for  which  a  suitable  boun¬ 
dary-conforming  coordinate  grid  can  be  construct¬ 
ed.  These  schemes  still  remain  quite  expensive 
in  terms  of  computer  resources  for  practical  use, 
however,  primarily  because  of  the  large  number  of 
grid  cells  necessary  for  adequate  resolution  in 
three-dimensional  problems  and  the  large  number 
of  iterations  required  to  achieve  even  modest 
convergence  on  these  fine  grids.  The  present 
paper  describes  work  addressed  primarily  at  this 
last  difficulty,  but  also  includes  an  improvement 
which  addresses  the  first  problem. 

The  major  thrust  of  the  current  work  is  the 
incorporation  of  a  multi-grid  algorithm4,5  to 
solve  the  difference  equations.  At  the  same 
time,  the  artificial  viscosity  terms  have  been 
modified  to  maintain  almost  everywhere  the 
second-order  accuracy  of  the  original  central- 
difference  approximation  used  in  subsonic  regions 
of  the  flowfield  in  a  manner  similar  to  that  used 
for  two-dimensional  calculations  by  Jameson  in 
Reference  6.  When  using  multi-grid  to  accelerate 
convergence  in  two-dimensional  calculations, 
Jameson7  found  it  necessary  to  use  a  generalised 
Alternating-Direction-Implicit  (ADI)  smoothing 
routine  to  eliminate  all  high  wavenumber  com¬ 
ponents  of  the  error,  however  the  results  of 
Shmilovich  and  Caughey8  and  their  extension  to 
the  current  work  demonstrates  that  good  rates  of 


convergence  can  be  obtained  using  modified  ver¬ 
sions  of  the  original  SLOR  algorithms.  In  order 
to  provide  reliable  convergence,  the  bandwidth  of 
the  original  SLOR  scheme  has  been  increased  to 
allow  pentadiagonal  inversions  along  each  line 
(instead  of  the  tridiagonal  inversions  of  the 
original  scheme). 

In  the  present  paper,  a  brief  review  of  the 
fully-conservative  finite  volume  scheme  will 
first  be  presented,  concentrating  upon  those 
aspects  necessary  for  an  understanding  of  the 
improvements  in  the  artificial  viscosity,  the 
modified  SLOR  schemes,  and  the  implementation  of 
the  multi-grid  algorithm.  The  changes  resulting 
from  the  implementation  of  the  new  features  will 
then  be  described,  and  results  indicating  the 
improved  iterative  convergence  and  accuracy  of 
the  new  scheme  will  be  presented.  A  comparison 
of  results  calculated  using  the  original  first- 
order  accurate  and  new  second-order  accurate 
schemes  will  also  provide  some  guidelines  on  the 
number  of  mesh  points  required  for  given  levels 
of  accuracy  in  force  coefficients  for  these 
three-dimensional  calculations. 


II.  ANALYSIS 

The  current  work  is  based  upon  the  finite- 
volume  method  of  Jameson  and  Caughey. 1-3  That 
method  provides  a  discrete  approximation  to  the 
nonlinear  potential  equation  of  transonic  flow 
which  may  be  interpreted  either  as  a  finite- 
difference  method  which  balances  fluxes  across 
cell  faces  or  as  a  finite  element  siethod  based 
upon  the  Bateman  variational  principle.  In  the 
original  formulation  of  that  method,  a  first- 
order  truncation  error  was  introduced  by  the 
addition  of  an  artificial  viscosity  needed  to 
stabilize  the  scheme  in  regions  of  supersonic 
flow,  and  the  difference  equations  were  solved  by 
an  SLOR  scheme. 


A.  Finite-Volume  Formulation 

Many  aerodynamic  problems  of  practical 
interest,  including  transonic  flows  with  weak 
snock  waves,  can  be  usefully  approximated  as 
potential  and  steady.  In  strong  conservation 
form,  the  equation  for  the  velocity  potential  • 
can  be  written  in  Cartesian  coordinates  (x,y,z) 
as 

(pVx  +  <«Vy  +  (PV,  ’  °>  (1) 
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where  p  is  the  density,  given  by  the  isentropic 

law 

p  -  {l  ♦  ( y-1)/2  -  q2)}1/(Y_1)  (2) 

Here  M  is  the  Mach  number  of  the  free  stream, 

q  is  the  magnitude  of  the  velocity  Vt,  and  the 
density  and  velocity  have  been  normalized  by 
their  freestreaa  values. 


Numbering  the  cell  vertices  as  illustrated 
in  Figure  1,  and  assuming  that  the  local 
coordinates  X.  ■  +1/2,  Y.  *  +  1/2,  Z.  *  +  1/2 

i  —  l  —  i  — 

at  the  vertices,  the  local  mapping  can  be  written 
4 

x  -  87  *.(1/4  +X.X)(l/4  +Y.Y)(l/4  +Z.Z).  (6) 

^ — 1  ii  i  i 

i-1 


The  finite-volume  method  for  transonic 
potential  flow1-3  is  a  geometrically-general 
technique  based  upon  a  numerical  evaluation  of 
the  transformation  aetrics  produced  by  an 
arbitrary  transformation  to  boundary-conforming 
coordinates.  Consider  a  transformation  to  a  new 
set  of  coordinates  X,Y,Z.  Let  the  Jacobian 
matrix  of  the  transformation  be  defined  by 


Similar  formulas  hold  for  y,  z  and  ♦  .  If  we 
introduce  the  averaging  and  differencing 
operators 


M„f.  .  .  “  l/2(f.  ...  .  +  f.  .  .  .) 

X  i,j,k  1+1/2, j,k  1-1/2, j,k 


«vf.  .  . 
^  * i J  »k 


(f i+1/2 , j ,k~  fi-l/2,j,k) 


then  the  transformation  dervivatives,  evaluated 
at  the  cell  centers,  can  be  expressed  by  formulas 
(3)  such  as 

x  “  u  6  x 
X  MYZ  X  ’ 


and  let  h  denote' the  determinant  of  H.  The 
metric  tensor  of  the  new  coordinate  system  is 

T 

given  by  the  matrix  G  ■  H  H  ,  and  the  con- 
travariant  components  of  the  velocity  vector 
U,  V,  and  W  are  given  by 


Eq.(l),  upon  multiplication  by  h,  can  then  be 
written 

(phU)x  ♦  (phV)y  ♦  (phW>2  -  0.  (5) 

The  fully-conservative,  finite-volume 
approximation  corresponding  to  Eq.(5)  is 
constructed  by  assuming  separate  trilinear 
variations  of  the  independent  and  dependent 
variables  within  each  mesh  cell. 


Z 


(a)  Physical  cell  (b)  Computational  cell 

Figure  1.  Mapping  of  mesh  cells. 


(8) 


with  similar  expressions  for  the  derivatives  of 
y,  z  and  the  potential.  Such  formulas  can  be 
used  to  determine  p,  h,  U,  V,  and  H  at  the 
center  of  each  cell  using  Eqs.(2),(3),  and  (4). 
Eq.(S)  is  represented  by  conserving  fluxes  across 
the  boundaries  of  auxiliary  cells  whose  faces  are 
chosen  to  be  midway  between  the  faces  of  the 
primary  mesh  cells.  This  can  be  represented  as 

U^SjiphU)  ♦  Wxz{Y<phV)  +  WXYSZ(phW>  “  °'  (9) 


This  formula  can  also  be  obtained  by  applying  the 
Bateman  variational  principle  that  the  integral 
of  the  pressure 

X  •  /  p  dx  dy  dz  (10) 


is  stationary,  and  approximating  I  by  a  simple 
one-point  integration  scheme  in  which  the 
pressure  at  the  center  of  each  grid  cell  is 
multiplied  by  the  cell  volume.  In  this  way,  for 
subsonic  flow,  the  finite-volume  method  can 
equally  well  be  regarded  as  a  finite  element 
method  with  isoparametric  trilinear  elements. 


The  use  of  the  one-point  integration  scheme 
leading  to  Eq.(9)  has  the  advantage  of  requiring 
only  one  density  evaluation  per  mesh  point,  but 
also  has  the  undesirable  effect  of  tending  to 
decouple  the  solution  at  odd-  and  even-numbered 
points  of  the  grid,  and  suitable  recoupling  terms 
can  be  added  to  improve  the  stability  of  the 
solut ion. 


We  define 


T  "  -fK5XY(AX  +  V^XY 


♦  Wx<yz<Ay  +  Ag>l*xfiYZ 
*  pysxz(ax  +  az>my*xz 

-1/26  (A  +A  +A  ) 6  }♦ 

XYZ  X  Y  Z  XYZ 
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where 


Ax  -  ph(g“  -  U2/a2), 

Ay  -  ph(g22  -  V2/a2),  (12) 

*z  “  ph(g33  -  W2/a2) , 


are  the  coefficients  of  A 


XX’ 


♦yy ,  and 

ij 


’ZZ’ 

are  the 
-1 


the  expanded  for*  of  Eq.(5).  Here  g 

elements  of  the  inverse  of  the  metric  tensor  G 
and  az  is  the  square  of  the  local  speed  of 
sound.  The  addition  of  T  to  Eq.(S)  provides 
recoupling  for  0  <  e  <  1/2.  For  e  *  1/2,  this 
reduces  Eq.(9)  to  the  usual  second-order 
accurate,  six-point  Laplacian  operator  for 
incompressible  flow  on  a  uniform  grid. 


P  -  v  P 
i.j.k  i,j,k  i-l.j.k* 

Pi+l/2,j,k"  ^  A 

-P .  ♦  v  P . 

l+l, j,k  l+l, J.k  i+2, j,k 


U  *  0, 
(18) 
U  <■  0, 


Similar  expressions  are  used  for  the  contributons 
from  the  Q  and  R  fluxes.  In  regions  where  the 
solution  is  smooth,  the  terra  ^x^YZ^  *8 

first  order  in  the  mesh  spacing,  and  the 
viscosity  is  formally  a  second  order  quantity. 
Near  a  shock,  for  an  appropriate  value  of  x,  the 

quantity  v.  .  becomes  small,  and  Eqs.(18) 
i  *  3  » k 

approximate  Eqs.(15)  —  i.e.,  the  viscosity 
reverts  to  the  original  first-order  form.  This 
hybridization  pf  the  second-order  scheme  has  been 
found  necessary  to  stabilize  computations  for 
solutions  containing  strong  shocks. 


The  scheme  is  stabilized  in  supersonic 
regions  by  the  explicit  addition  of  an  artificial 
viscosity.  The  viscosity  terms  added  in  the 
original  formulation  are  chosen  to  emulate  the 
directional  bias  introduced  by  the  rotated 
difference  scheme  of  JasKSon.  The  fluxes  $,  $, 
and  fi  are  defined  such  that 

Pi,j,k-  P^/«2(«2«xx  ♦  «V*xy  +  UW^)*, 

*i,j,k"  Pho/.2^  ♦  V2«yy  +  VW«yZ)*,  (13) 

*i,j,k“  Ph0/*2(UW8XZ  *  WSYZ  *  Az**’ 


where  the  switching  function 

a  -  max(0. ,  1  -  ( M  /M)2) 
c 


(14) 


is  non-zero  only  for  values  of  the  local  Mach 

number  M  greater  than  some  critical  Mach  number 

M  .  Then,  after  defining 
c 


p  a  iij>k 

*i+l/2,j,k  _  A 

l+l, j.k  ’ 


U  »  0, 

U  <  0, 


(15) 


with  similar  shifts  for  Q  and  R,  Eq.(9)  is 
represented  as 


«x(Myz(°hU)  ♦  P)  ♦  «y(uxz(phV)  ♦  Q) 
♦  VMXY(phW)  *  R)  ♦  T  -  0. 


(16) 


The  difference  Eqs.(9)  approximate  the 
original  differential  Eq.(5)  to  within  a  formal 
truncation  error  of  second  order  in  the  mesh 
spacing  in  the  physical  plane  when  the  mesh  is 
ssnoth.  Since  the  additional  fluxes  P,  Q,  and 
R  added  in  supercritical  regions  are  of  the 
order  of  the  physical  mesh  spacing  itself, 
however,  Eqs.(16)  approximate  Eq.(5)  to  within  a 
truncation  error  of  only  first  order  in  the  mesh 
spacing.  The  error  resulting  from  the  introduc¬ 
tion  of  the  artificial  viscosity  can  be  reduced 
to  second  order  nearly  everywhere  in  the  flow 
field  if  we  define* 


B.  Multi-Grid  Iteration 


The  difference  equations  resulting  from 
Eq.(16)  can  be  solved  by  carefully  constructed 
SLOR  schemes.  The  SLOR  schemes  described  in 
References  1-3  were  constructed  in  a  manner  that 
required  only  tridiagonal  inversions  along  each 
line.  When  the  contributions  arising  from  the 
inclusion  of  the  artificial  viscosity  terms  are 
included,  the  corrections  at  each  point  are  coup¬ 
led  to  those  of  its  two  neighbors  on  one  side 
(either  side  must  be  allowed  in  a  general  scheme, 
depending  u£on  the  sign  of  the  velocity)  for  the 
first-order  scheme,  and  itB  three  nearest  neigh¬ 
bors  on  one  side  for  the  second-order  form  of  the 
viscosity.  Thus,  a  general  scheme  which  accounts 
for  all  these  contributions  would  require  a  pen- 
tadiagonal  inversion  for  the  first-order  scheme, 
or  a  septadiagonal  inversion  for  the  second-order 
scheme.  It  was  found  that  the  pentadiagonal 
inversion  scheme  was  substantially  more  stable 
than  the  tridiagonal  scheme  when  the  second-order 
form  of  the  viscosity  was  used,  but  made  little 
difference  in  convergence  behavior  when  the 
first-order  viscosity  was  used.  Complementary 
experiments  by  A.  Jameson  of  Princeton  University 
have  shown  no  consistent  advantage  in  using  sep¬ 
tadiagonal  inversions  (over  the  pentadiagonal 
schesw)  when  the  second-order  viscosity  is  used. 
The  pentadiagonal  inversion  scheme  has  been 
incorporated  for  the  present  calculations. 

Both  X-line  and  Y-line  schemes  have  been 
implemented.  Only  the  Y-line  scheme  will  be 
described  here;  the  X-line  scheme  can  be  similar¬ 
ly  constructed.  Also,  the  coefficients  will  be 
described  only  for  the  case  when  U,V,W  >  0;  the 
coefficients  for  other  cases  can  easily  be 
constructed  by  analogy. 

We  define 

AB  -  w§PhU/a2  max(|u|,|v|,|w|), 

Ay-  u^phV/a2  max(|u|,|v|,|w|),  (19) 

Ay  ■  «gPhW/a2  max(|u|,|v|,|w|). 


Vi.j.k 


1  -  xJ^p 


(17) 


where  x  is  a  constant  of  order  unity,  and 


where  is  a  parameter  governing  the  amount  of 

•  type  damping  added  explicitly  to  the  time 


3 


*■  M 


§3 


*« 


dependent  equation  nodelling  the  relaxation 
proceaa.  Also, 

Ayu  -  Ph«2/«2, 

Ayy  -  phdV^a2, 

Ayy  -  phOW2/a2. 

Then  the  correction  to  the  potential 


Ci.  j.k 


■  C  n^  1 1 
t.  .  , 
i.J.k 


.(n) 

i.j.k 


is  calculated  according  to 

lAr*  'i.j-i,,' 

*  (AZ  +  j  ,k  ci,j,lt-l) 

+  (AX*  V<Ci.j,k’Ci-l,i.k) 


♦  A.  (C.  .  . 

*  1  >  j 


i^i.j y 


♦  \_J-C.  ,  .  .  ♦  (3  ♦  v.  .  )c.  .  .  (22) 

UUl  i*l,j,k  i.j.k  i,j,k 

-  (3  *2w.  .  .  )C.  ,  .  . 

i.J.k  i-l, J,k  l-l.j.k 

*  ^  *  Wi,j.k  *Wi-l,j,k)Ci-2,j,k^ 

♦  \vl(3  *Vi, j,k)Ci, j.k  (4  +3vi,j,k)Ci,j-l,k 

♦  \»l<3  *  ”i,J,k>ci,],k-il 


♦  (2/»  -  I)(AV  +  A  )C 

I  *•  l»Ji* 


.  . 

i.j.k 


where  R.  .  is  the  residual  of  Eq.(16), 
i.j.k 

calculated  uaing  values  of  the  potential  f row  the 
previous  iteration,  and  w  is  an  overrelaxation 
parameter,  which  is  set  to  2  in  supersonic 
regions.  Eq.(22)  requires  a  pentadiagonal 
inversion  along  each  i-line  since  for  U  <  0  the 
formula  wust  be  modified  to  include  the  effect  of 
the  correction  at  the  (i*2,j,k)  point.  The  V 
end  W  component a  should  be  non-negative  in 
superaonic  tones  for  the  relaxation  sweeps  to 
proceed  in  the  positive  T-  and  Z-  directions. 

Mote  that  the  influence  of  corrections  at  the 
(i,j*l,k)  and  (i,j,k+l),  (i,j,k-2)  points  as  well 
as  the  (i-3,j,k),  (i,j-3,k)  and  (i,j,k-3)  points 
have  been  eliminated  by  the  effective  addition  of 
■ixed  space- tine  differences. 

These  SLOR  schemes  have  the  advantages  of 
being  quite  stable,  and  of  rapidly  eliminating 
any  large  local  errors  in  the  initial  estimates 
for  the  potential  field.  Their  rates  of  conver¬ 
gence  decrease  as  the  local  errors  become  smal¬ 
ler,  however,  with  the.  result  that  convergence  to 
very  small  residuals  can  be  excruciatingly  slow, 
especially  when  the  swsh  spacing  is  small. 

An  efficient  alternative  has  been 
demonstrated  by  Jaamson,  based  upon  the 
multi-level  adaptive-grid  technique  first 
proposed  by  Fedorenko,  and  developed  and 
popularised  by  Brandt.  The  concept  behind  the 
multi-grid  method  is  to  eliminate  each  band  of 
wavenumbers  in  the  error  spectrus  on  a 


finite-difference  grid  which  is,  in  a  sense, 
optisuil  for  that  component.  Thus,  low  wavenumber 
errors  are  eliminated  on  coarse  grids,  and  only 
the  high  wavenumber  components  need  be  eliminated 
on  the  fine  grids.  Alternatively,  the  use  of 
coarse  grids  to  eliminate  the  low  wavenumber 
component  of  the  error  can  be  thought  of  as 
allowing  a  very  high  signal  speed  for  the  effects 
of  this  error  to  be  transmitted  across  the  grid. 

The  multi-grid  method  was  first  applied  to 
the  transonic  small  disturbance  equation  by  South 
and  Brandt,  who  noted  the  problems  associated 
with  highly  stretched  grids  when  using  SLOR  as 
the  smoothing  algorithm,  and  suggested  using 
alternating  SLOR  as  a  remedy.  Three-dimensional 
calculations  using  the  full  potential  equation 
(in  non-conservation  form)  have  been  performed  by 
McCarthy  and  Reyhner,1  for  the  transonic  flow 
past  axisymmetric  inlets.  Their  computations 
were  performed  in  a  non-body-aligned,  cylindrical 
coordinate  system. 

The  structure  of  the  multi-grid  method  is  as 
follows.  Let  the  discretization  of  Eq.(16)  be 
represented  by 


k  (n+1) 
L  v 


k-1,2 . K, 


on  a  heirarchy  of  grid  levels  G°,  G1 . G  , 

with  K  denoting  the  finest  grid.  The  iterative 
solution  is  started  from  some  initial  estimate  on 
the  finest  grid.  After  the  high  wavenumber 
component  of  the  error  has  been  eliminated,  the 
fine-grid  residual  is  calculated  and  restricted 
to  the  next  coarsest  grid.  On  all  but  the  finest 
grid,  the  residual  must  be  modified  to  account 
for  the  difference  in  truncation  error  on  the 

various  grid  levels  (i.e.,  L  ♦/  0  unless 
k  "  K,  where  •  is  the  converged  solution  on  the 
K-th  grid).  Thus 

pk-1  „  Lk-lf(n)  _  I|'-lLk#(n)>  (24) 

k-1 

where  1^  is  a  restriction  operator  which 

averages  the  residuals  over  the  fine  mesh  points 
in  the  vicinity  of  each  coarse  grid  point.  After 
the  high  wavenumber  error  on  the  coarser  grid  has 
been  eliminated,  the  finer  mesh  solution  can  be 
improved  according  to 


where  is  a  prolongation  operator  which  is 

used  to  interpolate  the  coarse  grid  corrections 
onto  the  finer  grid. 

While  the  essence  of  the  idea  has  been 
described  above  for  two  grid  levels,  the  idea  can 
be  extended  to  as  many  levels  as  feasible  in 
order  to  work  on  the  broadest  possible  band  of 
the  error  wavenumber  spectrum.  Useful  error 
reduction  can  be  achieved  on  very  coarse  grids, 
containing  only  a  few  cells  in  each  coordinate 
direction. 
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In  the  original  implementation  of  Brandt, 
an  adaptive  strategy  was  envisioned  for  deter¬ 
mining  when  to  change  from  one  grid  level  to  the 
next.  The  smoothing  would  continue  on  a  partic¬ 
ular  grid  level  until  the  convergence  rate  fell 
below  some  predetermined  tolerance,  at  which  time 
the  residual  from  that  grid  would  be  restricted 
to  the  next  coarsest  grid.  The  smoothing  on  that 
grid  would  proceed  until  the  convergence  rate 
again  slowed,  at  which  time  the  residual  would 
again  be  restricted  to  a  coarser  grid,  and  the 
process  repeated.  When  the  solution  had  conver¬ 
ged  on  the  coarsest  level,  the  corrections  would 
successively  be  added  back  to  the  finer  grid 
solutions,  and  the  cycle  would  be  repeated.  In 
the  present  implementation,  a  simple  fixed  stra¬ 
tegy  has  been  found  effective.  A  fixed  number  of 
relaxation  sweeps  is  performed  on  each  grid 
before  restricting  the  residual  from  that  grid  to 
the  next  coarsest  level,  and  a  fixed  number  of 
relaxation  sweeps  is  performed  on  each  grid  after 
the  corrections  are  added  from  the  preceding  grid 
before  the  corrections  are  added  to  the  next 
finest  grid. 

In  the  present  codes,  as  in  our  earlier 
work,8  the  restriction  operator  averages  the 
residuals  at  the  27  fine  grid  neighbors  of  each 
coarse  grid  point,  weighted  according  to  the 
fraction  of  the  coarse  grid  cell  volume  asociated 
with  each  fine  grid  cell  in  the  computational 
space.  The  prolongation  operator  uses  four-point 
Lagrangian  interpolation  in  each  coordinate 
direction,  except  near  boundaries  where  the  order 
must  be  reduced. 

The  computational  labor  required  for  one 
multi-grid  cycle  using  this  fixed  strategy  is  as 
follows.  Let  one  work  unit  be  defined  as  the 
labor  required  for  one  relaxation  sweep  on  the 
finest  grid.  Then  if  mt  sweeps  are  done  on 
each  grid  before  the  residual  is  restricted  to 
the  next  coarsest  grid,  and  m2  sweeps  are  done 
on  each  grid  before  the  corrections  are  prolonged 
to  the  next  finest  grid,  then  the  cost  of  one 
cosiplete  multi-grid  cycle  is  approximately 

{Bim,  ♦  1)  +  m2}/7  work  units. 

This  includes  the  cost  of  computing  the  residual 
on  each  grid  for  restriction  to  the  next  coarsest 
grid  as  approximately  equal  that  of  a  relaxation 
sweep  on  that  grid  (since  most  of  the  labor  is 
involved  in  computing  the  residual,  not  in  the 
actual  update  of  the  solution),  but  neglects  the 
overhead  in  restricting  residuals  and  prolonging 
corrections.  For  m4  ■  mg  *  1,  one  multi-grid 
cycle  requires  approximately  2-3/7  work  units. 

The  success  of  the  multiple-grid  method 
depends  upon  the  efficient  elimination  of  high 
wavenumber  errors  on  any  given  grid.  Jameson8 
used  a  generalized  alternating-direction  scheme, 
in  which  he  replaced  the  usual  constant  in  each 
factor  by  the  sum  of  a  constant  and  first-order 
difference  operators  in  each  coordinate  direc¬ 
tion.  Shmilovich  and  Caughey8  have  shown  that, 
even  for  SLOR  schemes,  the  growth  factor  in  a  von 
Neumann  analysis  should  never  exceed  approxi¬ 
mately  0.78  per  multi-grid  cycle  if  the  multi¬ 
grid  algorithm  is  effective  on  error  with  low 
wavenumber  components  in  any  of  the  three 


coordinate  directions.  The  effectiveness  of  the 
multi-grid  procedure  in  eliminating  error  having 
low  wavenumber  component  in  only  one  (or  two) 
directions  was  not  verified  in  the  present  work, 
but  it  was  found  effective  to  alternate  between 
X-line  and  Y-line  SLOR  in  conjunction  with  the 
multi-grid  procedure.  This  can  be  done  in  two 
different  ways:  (1)  alternate  multi-grid  cycles 
can  be  performed  using  the  two  schemes,  or 
(2)  the  two  schemes  can  be  alternated  at  each 
level  within  each  multi-grid  cycle  (if  nij  and 
m2  are  greater  than  one).  The  most  effective 
procedure  seems  to  be  the  second  option  (with 
m4  ■  2  and  ra2  «  4  or  6). 

C.  Geometrical  Aspects 

The  algorithm  described  above  has  been 
incorporated  into  two  computer  programs  for 
calculating  the  transonic  potential  flow  past 
three-dimensional  wings  and  wing-fuselage  combi¬ 
nations.  These  codes  are  known  generally  as 
FLO-27  and  FLO-30.  FLO-27  analyses  the  flow 
past  swept  wings  of  essentially  arbitrary  plan- 
form  and  section  shape;  FLO-30  analyses  the 
flow  past  such  wings  mounted  upon  arbitrary 
fuselage  shapes. 

Both  codes  construct  boundary-conforming 
coordinate  grids  by  sequences  of  simple  conformal 
and  shearing  transformations.  The  computational 
domain  in  each  program  is  terminated  at  artifi¬ 
cial  boundaries,  located  approximately  ten  chords 
distant  from  the  wing  surface  in  each  spanwise 
surface,  and  approximately  four  semi-spans  from 
the  symmetry  plane  or  fuselage  in  the  lateral 
direction.  On  the  upstream  and  lateral  boun¬ 
daries,  the  potential  describing  perturbations 
from  the  uniform  free  stream  is  set  to  zero, 
while  on  the  downstream  boundary,  the  velocity 
perturbations  in  the  streamwise  direction  are  set 
to  zero  (consistent  with  a  fully-developed  flow 
in  the  Trefftz  plane).  The  no-flux  condition  is 
enforced  directly  in  the  flux  balances  at  solid 
boundaries,  and  a  linearized  approximation  to  the 
vortex  sheet,  which  assumes  the  shed  vorticity 
trails  in  the  freestream  direction  in  a  fixed 
surface  downstream  of  the  trailing  edge  is  used. 
The  flux  balance  represented  by  Eq.(16)  is  also 
satisfied  at  points  on  the  vortex  sheet,  since  it 
does  not  require  differences  across  the  sheet. 

These  codes,  and  their  associated  grid 
generation  techniques,  are  described  in  greater 
detail  in  References  1-3. 


III.  RESULTS 

Results  will  be  presented  illustrating  both 
the  improved  rate  of  convergence  of  the  multi¬ 
grid  algorithm,  and  the  increased  accuracy  of  the 
scheme  with  the  second-order  viscosity. 

A.  Computational  Aspects 

Both  programs  have  been  designed  to  run  on 
either  modest  computers  with  large  virtual  memory 
or  on  advanced  machines  with  large  high-speed 
memories.  Even  so,  only  the  Cartesian 
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coordinates  of  the  mesh  and  the  solution  vector 
can  reasonably  be  stored  for  fine  meshes,  and  the 
transformation  derivatives  are  recomputed  at  each 
mesh  point  in  each  iteration.  Largely  because  of 
this,  the  program  runs  efficiently  on  vector 
machines  even  though  the  actual  line  inversions 
for  the  solution  are  inherently  recursive.  On  a 
grid  containing  160x24x32  mesh  cells  in  the 
X,  Y,  and  Z  directions,  respectively,  FLO- 30 
requires  about  830,000  words  of  storage  on  the 
Cray-I  computer,  and  requires  about  2.0  CPU 
seconds  per  work  unit  on  this  grid.  This  corres¬ 
ponds  to  an  estimated  average  computational  rate 
of  32  megaflops.  A  typical  solution  is  conver¬ 
ged  to  within  reasonable  engineering  accuracy 
after  about  30  work  units;  this  requires 
approximately  63  seconds  on  this  (relatively 
fine)  grid. 


B.  Computed  Results 

The  first  results  to  be  presented  are  for 
the  high-aspect-ratio,  supercritical  wing  (Wing 
A)  tested  by  Hinson  and  Burdges.11  A  perspective 
view  of  the  wing  is  shorn  in  Figure  2.  The  flow 
past  the  wing  at  a  free  stream  Mach  number  of  0.82 
and  1.3  degrees  angle  of  attack  was  analysed 
using  FLO-27.  The  analysis  was  performed  on  a 
sequence  of  three  grids,  each  obtained  by  doub¬ 
ling  the  number  of  mesh  cells  in  each  coordinate 
direction  from  the  preceding  grid,  and  prolonging 
the  results  of  the  converged  solution  from  the 
preceding  grid  as  the  initial  estisiate  on  the 
next  grid.  The  finest  grid  contained  128x16x32 
mesh  cells.  Calculations  were  performed  with 
both  the  first-  and  second-order  forms  of  the 
artificial  viscosity;  the  iterative  convergence 
rates  were  nearly  identical.  Figure  3  shows  the 
convergence  history  of  the  second-order  scheme  on 
the  finest  grid.  The  logarithm  of  the  average 
residual,  the  root  section  lift  coefficient,  and 
the  number  of  supersonic  points  are  plotted  as  a 
function  of  computational  work  (measured  in  work 
units);  the  lift  coefficient  is  normalized  by 
its  final  converged  value,  and  number  of  super¬ 
sonic  points  is  normalized  by  twice  its  final 
converged  value,  while  the  residual  is  normalized 
by  its  initial  value.  The  solid  lines  represent 
the  convergence  of  the  multi-grid  algorithm  (with 
mi  ■  2,  mj  "  6,  and  the  alternating  SLOR  scheme) 
using  4  grid  levels,  and  the  dashed  lines  rep¬ 
resent  the  convergence  of  a  relaxation  solution 
for  the  same  initial  quess.  Note  that  after  even 
100  relaxation  sweeps,  the  SLOR  scheme  has  elim¬ 
inated  only  about  half  of  the  error  in  root  sec¬ 
tion  lift  coefficient  and  in  the  number  of  super¬ 
sonic  points.  This  illustrates  the  slowness  with 
which  SLOR  eliminates  the  low  wavenumber  com¬ 
ponent  of  error.  With  the  multi-grid  scheme, 
both  of  these  measures  have  converged  to  within 
the  plottable  accuracy  of  the  figure  in  the 
equivalent  of  50  relaxation  sweeps.  The  wing 
surface  pressure  distribution  for  the  first-  and 
second-order  accurate  schemes  are  presented  in 
Figures  4(a)  and  4(b),  and  the  streamwise  pres¬ 
sure  distributions  at  the  25  per  cent  semi-span 
station  are  presented  for  both  schemes  in  Figures 
3(a)  and  5(b).  Note  the  increased  sharpness  with 
which  the  shocks  are  resolved  by  the  higher-order 
scheme.  Finally,  a  convergence  study  of  the  wing 
lift  and  drag  coefficients  with  mesh  specing  is 


shown  in  Figure  6.  Both  first-  and  second-order 
accurate  results  are  plotted;  the  former  vs. 
mesh  spacing  and  the  latter  ^s.  the  square  of  the 
mesh  spacing.  Straight  lines  through  the  finer 
mesh  results  for  both  schemes  converge  to  the 
same  asymptotic  value  for  the  lift  coefficient  in 
the  limit  of  zero  mesh  width,  but  the  drag 
results  have  not  yet  reached  their  asymptotic 
rates  on  these  grids.  The  absolute  error  in  lift 
for  the  second-order  scheme  on  the  finest  grid  is 
about  2  per  cent,  while  a  meoh  containing  more 
than  eight  times  as  many  cells  (a  factor  of  2  in 
each  coordinate  direction)  would  be  required  for 
similar  accuracy  using  the  first-order  scheme. 

The  absolute  error  in  the  drag  coefficient  for 
the  second-order  scheme  is  about  3  per  cent , 
while  similar  accuracy  for  using  the  first-order 
scheme  would  seem  to  require  approximately  64 
times  as  many  mesh  cells  (a  factor  of  four  in 
each  coordinate  direction). 

Finally,  to  illustrate  the  reliability  of 
the  scheme,  results  for  a  strongly  supercritical 
case  a^e  presented.  The  flow  past  the  ONERA  M-6 
wing,1  mounted  upon  a  circular  cylinder,  was 
computed  using  FLO-30.  Figure  7  shows  the  coor¬ 
dinate  lines  in  the  wing  and  fuselage  surfaces 
for  the  grid  used;  only  every  fourth  line  is 
shown  for  clarity.  The  solution  was  computed  at 
a  freestream  Mach  number  of  0.923  and  3.06 
degrees  angle  of  attack  on  a  very  fine  grid  con¬ 
taining  160x24x32  mesh  cells.  Nearly  20  per 
cent  of  the  mesh  points  have  supersonic  veloc¬ 
ities  for  this  case.  The  wing  pressure  distribu¬ 
tion,  showing  the  strong  shocks  at  the  trailing 
edge  of  the  upper  surface  and  the  substantial 
supersonic  pocket  outboard  on  the  lower  surface, 
is  plotted  in  Figure  8.  The  convergence  history 
is  plotted  in  Figure  9.  Again,  the  root  section 
lift  coefficient  and  the  number  of  supersonic 
points  have  converged  to  within  plottable  accu¬ 
racy  in  the  equivalent  of  about  50  relaxation 
sweeps . 


IV.  CONCLUDING  REMARKS 

A  multi-grid  algorithm  has  been  combined 
with  a  successive-line-overrelaxation  (SLOR) 
iterative  scheme  to  provide  improved  rates  of 
convergence  in  the  iterative  sense  for  the  com¬ 
putation  of  transonic  potential  flows  past  swept 
wing  and  wing-fuselage  configurations.  At  the 
same  time,  a  modified  form  of  artificial  viscos¬ 
ity  has  been  incorporated  which  results  in 
second-order  accuracy  for  the  scheme  nearly 
everywhere  in  the  flow  field.  The  method  has 
been  incorporated  into  two  computer  programs  for 
calculating  the  transonic  potential  flow  past 
three-dimensional  wings  and  wing-fuselage  com¬ 
binations.  Results  indicate  that  convergence 
adequate  for  most  engineering  purposes  can  be 
achieved  with  the  new  multi-grid  algorithm  in 
less  than  the  time  required  for  about  50  relax¬ 
ation  sweeps  using  the  original  SLOR  scheme. 
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Figure  2.  Geometry  of  Lockheed  Wing  A;  sections 
at  computational  stations  on  fine  grid  are  shown. 
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Figure  3.  Iterative  convergence  for  calculation 
of  flow  past  Lockheed  Wing  A  at  ■  0.82  and 
1.5  degrees  angle  of  attack;  second-order 
scheme . 
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Figure  4.  Three-dimensional  wing  surface  pressure  distributions  for  Lockheed 
Wing  A  at  M  “  0.82  and  1-5  degrees  angle  of  attack. 
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Figure  5.  Streaowise  surface  pressure  distributions  for  Lockheed  Wing  A  at  25 
per  cent  seai-span  station.  Sane  freestream  conditions  as  Figure  4. 


Figure  6.  Convergence  study  of  wing  lift  and 
drag  coefficients  for  Lockheed  Wing  A.  Same 
freestream  conditions  as  Figures  4  and  5. 

Figure  7.  ONERA  wing  M-6  mounted  upon  cylin¬ 
drical  fuselage.  Grid  lines  shown  in  wing  and 
fuselage  surfaces  only;  every  fourth  line  shown. 
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Figure  8.  Wing  surface  pressure  distribution  for 
flow  past  ONERA  Wing-cylinder  combination  at 
M"  "  0.923  and  3.06  degrees  angle  of  attack. 


Figure  9.  Iterative  convergence  for  calculation 
of  flow  past  ONERA  wing-cylinder  combination. 
Freestream  conditions  as  in  Figure  8. 
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